Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Speed up orbits of permutation groups on integers #4307

Open
wants to merge 5 commits into
base: master
Choose a base branch
from

Conversation

mjrodgers
Copy link
Collaborator

This aims to address issue #3287 by letting GAP know when we are in certain special cases which have optimized GAP methods for computing orbits.

Computing the orbit of a single point seems to still be slower than calling GAP directly by about a factor of 2, but this is a big improvement to the previous case (where it was about a 6x slowdown).

@mjrodgers
Copy link
Collaborator Author

@ThomasBreuer If you can take a look. At the moment I am only mapping ^ to OnPoints, I have some of the other conversions commented out since I don't think they are things that are supported yet.

src/Groups/gsets.jl Outdated Show resolved Hide resolved
src/Groups/gsets.jl Outdated Show resolved Hide resolved
@fingolfin
Copy link
Member

What about other places passing an action function to GAP?

@mjrodgers
Copy link
Collaborator Author

What about other places passing an action function to GAP?

I don't think there are any, when I checked. We don't currently experience any slowdown like this when calling stabilizer, for example. (I did a project-wide search for GapObj(action_function and this was the only place it was used.)

@ThomasBreuer
Copy link
Member

I think we have to distinguish two situations:

  • If the elements of Omega can be translated to pure GAP objects and if we have a Julia action function such that the GapObjs of the group generators act on the GAP objects in question via this function then we have a chance to hit an efficient GAP method for the computations, and we can first translate the input to GAP and in the end carry back the result to Julia.
    This is the case for permutation groups acting on positive integers or vectors or sets of them.
    (For the action on sets and vectors of integers, we have to convert also omega to a pure GAP object.)
  • In general we have some group acting on some objects via some function, for example a permutation group of degree n acting on polynomials in n variables by permuting the variables.
    (For examples, see test/Groups/gsets.jl).
    Here the old implementation uses that GAP allows us to act with Oscar group elements on Oscar objects, via the Julia function wrapped into a GAP function.

@ThomasBreuer
Copy link
Member

Do we have a better idea than keeping the generic orbit method and adding methods for those special cases that "can be translated to GAP and back", similar to the approach that was taken for stabilizer in #4281?

@fingolfin
Copy link
Member

Tests fail.

Also what this really is missing is some systematic benchmark tests, so that we can verify each change for how it affects performance.

Also, could you please share your example that where this PR makes it go from a 6x to a 2x slow down?

@mjrodgers
Copy link
Collaborator Author

I've added a new test to specifically catch the error (which was previously thrown outside of a @test statement). This case seems to throw an error when we use acts = GapObj(gens(G); recursive=true), but not using acts = GapObj(gens(G)); on the other hand, we get an error computing the orbits of symmetric_group(5) unless we use recursive=true... @ThomasBreuer any easy fix for this?

@mjrodgers
Copy link
Collaborator Author

mjrodgers commented Nov 20, 2024

Working to put together some systematic benchmarks, but for this specific case we have

julia> G = symmetric_group(5)
Sym(5)

julia> @benchmark Vector{Int}(GAP.Globals.Orbit(G.X,1))
BenchmarkTools.Trial: 10000 samples with 10 evaluations.
 Range (min … max):  1.700 μs …  4.329 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.829 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.844 μs ± 86.960 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

            ▁▄▆█▇▅▆
  ▂▂▂▂▂▃▃▄▇▇███████▇▆▅▅▅▄▄▅▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂ ▃
  1.7 μs         Histogram: frequency by time        2.21 μs <

 Memory estimate: 1.09 KiB, allocs estimate: 28.

With this pull request:

julia> @benchmark collect(orbit(G,1))
BenchmarkTools.Trial: 10000 samples with 8 evaluations.
 Range (min … max):  3.688 μs …   8.714 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     4.146 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   4.167 μs ± 195.265 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

                   ▄▇▄█▇▆▆▆▇▃▂
  ▁▁▁▂▂▂▃▃▂▂▂▂▂▂▃▅█████████████▅▇▅▄▃▃▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  3.69 μs         Histogram: frequency by time        4.89 μs <

 Memory estimate: 3.92 KiB, allocs estimate: 91.

This 2x slowdown is consistent when using larger symmetric_groups as well.

@@ -351,8 +362,8 @@ julia> length(orbit(Omega, 1))
"""
function orbit(Omega::GSetByElements{<:GAPGroup, S}, omega::S) where S
G = acting_group(Omega)
acts = GapObj(gens(G))
gfun = GapObj(action_function(Omega))
acts = GapObj(gens(G); recursive=true)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
acts = GapObj(gens(G); recursive=true)
acts = [GapObj(g) for in gens(G)]

but even this will fail if g cannot be converted to a GAP object...

@ThomasBreuer
Copy link
Member

ThomasBreuer commented Nov 20, 2024

Sorry if my comment above ("I think we have to distinguish two situations.") was not clear enough.

What I wanted to say is that in the current setup, we need first of all one generic orbit function. Up to now, we have one that delegates to GAP but keeps the Julia objects, that is, there is a GAP list of Julia group elements acting with a Julia function on Julia objects.
This method can deal with all cases but is expected to be slow.
(For this generic situation, an orbit method written in pure Julia is expected to be faster.)

In addition to that method, we can install more orbit methods. They can delegate to GAP or be written in Julia.
In the former case, they can decide to translate their arguments entirely to GAP, and then to translate the result back to Julia; this applies in particular to the situations for which GAP has special methods, such as the action of permutations on integers or (nested) lists or sets of them. (In pull request #4281, this has been done for stabilizer, it would look similar for orbit.) It is still not clear whether also in these situations, an orbit written in Julia will be faster than delegating to GAP.

@fingolfin
Copy link
Member

This still fails all its tests. @mjrodgers what's the status there, are you working on it?

Regarding the benchmarks: the first example actually suffers from inefficient code for converting a GAP list of integers into a Vector{Int}, so the discrepancy is actually even larger if one considers that this conversion currently is very suboptimal (see oscar-system/GAP.jl#777 for stuck WIP on improving this -- perhaps @ThomasBreuer will adopt this at some point? he already took care of the other conversion direction)

Also type instability hurts this micro benchmark. To wit:

julia> f(obj::GapObj) = Int[obj[i] for i in 1:length(obj)]  # better conversion code
f (generic function with 1 method)

julia> G = symmetric_group(5)
Sym(5)

julia> @benchmark Vector{Int}(GAP.Globals.Orbit(G.X,1))
BenchmarkTools.Trial: 10000 samples with 9 evaluations.
 Range (min … max):  2.361 μs …  15.579 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.440 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.463 μs ± 178.038 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

     ▂▆▇▇█▆▅▄▃▁
  ▁▂▆███████████▇▅▅▄▃▃▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  2.36 μs         Histogram: frequency by time         2.9 μs <

 Memory estimate: 1.09 KiB, allocs estimate: 28.

julia> @benchmark Vector{Int}(GAP.Globals.Orbit(G.X,1)::GapObj)
BenchmarkTools.Trial: 10000 samples with 9 evaluations.
 Range (min … max):  2.199 μs …  12.412 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.296 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.315 μs ± 164.926 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

         ▃▆█▇▅▂
  ▁▂▄▇▇▇███████▇▅▅▄▃▃▃▃▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  2.2 μs          Histogram: frequency by time        2.75 μs <

 Memory estimate: 1.09 KiB, allocs estimate: 28.

julia> @benchmark f(GAP.Globals.Orbit(G.X,1))
BenchmarkTools.Trial: 10000 samples with 9 evaluations.
 Range (min … max):  2.130 μs …  12.282 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.204 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.220 μs ± 149.329 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

         ▂▆▁█ ▄ ▁
  ▂▃▃▇▅▇▅██████▆█▄▅▄▃▃▂▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▁▂▂▂▂▂▂ ▃
  2.13 μs         Histogram: frequency by time        2.56 μs <

 Memory estimate: 720 bytes, allocs estimate: 21.

julia> @benchmark f(GAP.Globals.Orbit(G.X,1)::GapObj)
BenchmarkTools.Trial: 10000 samples with 9 evaluations.
 Range (min … max):  2.116 μs …  14.546 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.199 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.215 μs ± 181.758 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

              █▁▁▁▆
  ▁▁▂▂▅▃▄▄█▇▇██████▆▅▅▄▆▂▂▂▃▂▂▂▃▂▂▂▂▂▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  2.12 μs         Histogram: frequency by time        2.45 μs <

 Memory estimate: 720 bytes, allocs estimate: 21.

Of course this is minor and it equally affects the Oscar wrapper for GAP's Orbit.

For that, as @ThomasBreuer already pointed out, and as we discussed a few days ago, you'll probably need a different strategy:

E.g. it only makes sense to use the "the underlying GAP objects of the group generators" if the given group actually is a GAP group (which the method signature already ensures), and the objects being acted on are "native" GAP objects -- which right now probably means only for Int.

To unwrap the generators it is probably more efficient to use GAPWrap.GeneratorsOfGroup(GapObj(G)). Indeed (note the use of $G instead of G to get accurate results):

julia> G = symmetric_group(5);

julia> @benchmark GapObj(gens($G); recursive=true)
BenchmarkTools.Trial: 10000 samples with 209 evaluations.
 Range (min … max):  363.239 ns …  27.510 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     400.718 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   409.789 ns ± 284.250 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

            ▄▅█▆▅▅▂▂▂
  ▂▃▅▅▆▄▄▄▅███████████▇█▇▇▇▆▆▅▄▃▃▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  363 ns           Histogram: frequency by time          513 ns <

 Memory estimate: 192 bytes, allocs estimate: 6.

julia> @benchmark Oscar.GAPWrap.GeneratorsOfGroup(GapObj($G))
BenchmarkTools.Trial: 10000 samples with 998 evaluations.
 Range (min … max):  14.946 ns … 35.363 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     15.030 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   15.167 ns ±  0.742 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▃▇█▅▅ ▄▂ ▂ ▃ ▁▁▂▁                                           ▂
  █████▁████▁█████████▇▁▇███▁▆▇▆▆▆▁▇▆▅▆▁▆▅▅▆▁▅▄▄▅▄▁▁▃▅▅▁▄▄▅▅▄ █
  14.9 ns      Histogram: log(frequency) by time        17 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

@fingolfin
Copy link
Member

Next I wonder how expensive storing the orbit as an attribute is -- this potentially allocates a fresh Dict which can be costly. You could find out by temporarily disabling the set_attribute! and measuring again

@mjrodgers
Copy link
Collaborator Author

This still fails all its tests. @mjrodgers what's the status there, are you working on it?

Working with Thomas on it, this pull request might end up being superseded by the solution that Thomas describes in his previous comment; but I wanted to at least include this additional test here, since it can help track an extra "weird" case that might need some extra care.

@ThomasBreuer
Copy link
Member

Shall I make a new pull request that adds some special orbit methods, as I had sketched above, or shall I commit to the current pull request?

@mjrodgers
Copy link
Collaborator Author

I think it is good to add to the current pull request, if that is easy.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants